---
title: "Replication script for the article Experience with Plant-Based Meat Substitutes Promotes Sustainable Diets and the Intention to Eat Less Meat"
author: "Maiken Maier, Fabian Baettig, Leon Sistek, and Ursa Bernardic"
date: "18-03-2026"
output: 
  html_document:
    theme: paper #specifies the theme (structure/front)
    highlight: zenburn #creates dark background for code
    df_print: paged
    toc: True #table of contents
---

```{r prep_data, warning=FALSE,echo=FALSE, results=FALSE, message=FALSE}
library(readxl)
library(tidyverse)
library(writexl)
library(dplyr)
library(naniar)
library(psych)
library(jtools)
library(sandwich)
library(lmtest)
library(ggstance)
library(broom.mixed)
library(ggpubr)
library(effsize)
library(stats) #contains p.adjust formula
library(stringr)
library(openxlsx)
library(ggplot2)
library(gridExtra)
library(grid)
library(knitr)
library(kableExtra)
library(MASS) # For polr function
library(readr)
library(lavaan) #mediation model package
library(stargazer)
library(broom)
library(mediation)
library(estimatr)
library(modelsummary)
library(patchwork)
library(effectsize)
library(interactions)
library(grid)
library(estimatr)
library(modelsummary)
library(WeightIt)
library(cobalt)
```

# Data Cleaning & Recoding for Analysis
```{r}
## load data
tot_df <- read.csv('/Users/maimaier/Library/CloudStorage/OneDrive-ETHZurich/Conferences and Papers/Maier et al. 2025_Substitute survey paper/Submission/Replication data and code/SubstituteStudyI_replication_data.csv')
table(tot_df$PG_Chicken_quantity)
```


## Recode/clean merged dataset
```{r}
## Option 1: Generate product treated variable to explore subgroup effects
tot_df <- tot_df %>% mutate(product_treated = case_when(group == 0 ~ "control",
                                                        group == 1 & guggeli == 1 & kebab == 0 & skewer == 0 ~ "guggeli_only",
                                                        group == 1 & guggeli == 0 & kebab == 1 & skewer == 0 ~ "kebab_only",
                                                        group == 1 & guggeli == 0 & kebab == 0 & skewer == 1 ~ "skewer_only",
                                                        group == 1 & guggeli == 1 & kebab == 1 & skewer == 0 ~ "guggeliANDkebab",
                                                        group == 1 & guggeli == 1 & kebab == 0 & skewer == 1 ~ "guggeliANDskewer",
                                                        group == 1 & guggeli == 0 & kebab == 1 & skewer == 1 ~ "kebabANDskewer",
                                                        group == 1 & guggeli == 1 & kebab == 1 & skewer == 1 ~ "all_products"
                                                        ))
table(tot_df$product_treated)

## Option 2: Generate data subsets for individuals who sampled specific products to explore subgroup effects
skewer_df <- tot_df %>%
  filter(group == 0 | skewer == 1) %>%
  mutate(group = ifelse(skewer == 1, "skewer", "control"))
table(skewer_df$group)

guggeli_df <- tot_df %>%
  filter(group == 0 | guggeli == 1) %>%
  mutate(group = ifelse(guggeli == 1, "guggeli", "control"))
table(guggeli_df$group)

kebab_df <- tot_df %>%
  filter(group == 0 | kebab == 1) %>%
  mutate(group = ifelse(kebab == 1, "kebab", "control"))
table(kebab_df$group)

## Generate factor variables for mediation analyses
# Meat consumption
table(tot_df$meat_consumption)
tot_df <- tot_df %>% mutate(meat_consumption_f = case_when(meat_consumption <= 3 ~ "low_meat",
                                                           meat_consumption == 4 ~ "moderate_meat",
                                                           meat_consumption >= 5 ~ "high_meat"
                                                           ))
table(tot_df$meat_consumption_f)
tot_df$meat_consumption_f <- factor(tot_df$meat_consumption_f) #convert to factor
tot_df$meat_consumption_f <- relevel(tot_df$meat_consumption_f, ref = "low_meat") #Set baseline level

# Substitute consumption
table(tot_df$substitute_consum)
tot_df <- tot_df %>% mutate(substitute_consum_f = case_when(substitute_consum == 0 ~ "no_sub",
                                                            substitute_consum == 1 | substitute_consum == 2 ~ "low_sub",
                                                            substitute_consum >= 3 | substitute_consum == 4 ~ "moderatetohigh_sub"
                                                            ))
table(tot_df$substitute_consum_f)
tot_df$substitute_consum_f <- factor(tot_df$substitute_consum_f) #convert to factor
tot_df$substitute_consum_f <- relevel(tot_df$substitute_consum_f, ref = "no_sub") #Set baseline level

# Substitute opinion WI
table(tot_df$sub_general_opinion_w1)
tot_df <- tot_df %>% mutate(sub_general_opinion_w1_f = case_when(sub_general_opinion_w1 <= 3 ~ "neg_opinion",
                                                                 sub_general_opinion_w1 == 4 ~ "neutral_opinion",
                                                                 sub_general_opinion_w1 >= 5 ~ "pos_opinion"
                                                                 ))
table(tot_df$sub_general_opinion_w1_f)
tot_df$sub_general_opinion_w1_f <- factor(tot_df$sub_general_opinion_w1_f) #convert to factor
tot_df$sub_general_opinion_w1_f <- relevel(tot_df$sub_general_opinion_w1_f, ref = "neg_opinion") #Set baseline level
```


# Main Analysis with Control Variables
## Main Regression Models with Controls
```{r}
# Regressionen "group" mit Controls
dv_control_quality_index <- lm_robust(quality_index ~ group + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_health_index <- lm_robust(health_index ~ group + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_convenience_index <- lm_robust(convenience_index ~ group + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_environment_index <- lm_robust(environment_index ~ group + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_animal_index <- lm_robust(animal_index ~ group + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_price_index <- lm_robust(price_index ~ group + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_policy_support_1 <- lm_robust(policy_support_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_2 <- lm_robust(policy_support_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_3 <- lm_robust(policy_support_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)


dv_control_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

# Export regression table to LaTeX - Intention Variables
modelsummary(list("DV2a: Reduce meat consumption in the future" = dv_control_meat_intentions_1,
                  "DV2b: Reduce meat consumption in the next two weeks" = dv_control_meat_intentions_2,
                  "DV3a: Purchase PBMS in the future" = dv_control_sub_intentions_1,
                  "DV3a: Purchase PBMS in the next two weeks" = dv_control_sub_intentions_2),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results II (intention variables)",
  escape = FALSE
)

# Export regression table to LaTeX - Policy Support Variables
modelsummary(list("DV5a: Policies expanding plant-based food production/consumption" = dv_control_policy_support_1,
                  "DV5b: Policies promoting PBMS" = dv_control_policy_support_2,
                  "DV5c: Policies promoting meat reduction in CH" = dv_control_policy_support_3),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results III (general policy support)",
  escape = FALSE
)

# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1a: Quality" = dv_control_quality_index,
                  "DV1b: Health" = dv_control_health_index,
                  "DV1c: Convenience" = dv_control_convenience_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results I (PBMS evaluation indexes)",
  escape = FALSE
)
  
# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1d: Environment and Climate" = dv_control_environment_index,
                  "DV1e: Animal welfare" = dv_control_animal_index,
                  "DV1f: Price" =dv_control_price_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results I (PBMS evaluation indexes - continued)",
  escape = FALSE
)
```

## Stated-choice task: Product Shares

```{r}
# Helper: use rlang's %||% if available, otherwise define it
if (!exists("%||%")) {
  `%||%` <- function(x, y) if (is.null(x)) y else x
}

# ------------------------------------------------------------
# PRICE LISTS (per package)
# ------------------------------------------------------------
if (!exists("prices")) {
  prices <- list(
    Migros = list(
      PG_Chicken = 9.40, QP_Chicken = 10.55, NP_Chicken = 11.45,
      PG_Beef = 11.20, QP_Beef = 16.00, NP_Beef = 13.50,
      Planted_Kebab = 5.95, Planted_Guggeli = 5.95, Planted_Skewers = 6.95,
      Tofu_PG = 1.55, Tofu_Karma = 3.30, Tofu_Bio = 3.50
    ),
    Coop = list(
      PG_Chicken = 7.05, QP_Chicken = 12.00, NP_Chicken = 12.30,
      PG_Beef = 11.70, QP_Beef = 15.80, NP_Beef = 13.10,
      PG_ChickenMarinated = 2.95, Bell_ChickenMarinated = 7.90, NP_ChickenMarinated = 11.10,
      GM_Chicken = 5.95, GM_Beef = 5.95, GM_ChickenMarinated = 5.95,
      Planted_Kebab = 5.95, Planted_Guggeli = 5.95, Planted_Skewers = 6.95,
      Tofu_PG = 1.60, Tofu_Karma = 2.95, Tofu_Bio = 3.95
    )
  )
}

# ------------------------------------------------------------
# 1) CALCULATE STATED-CHOICE SHARES (cost-weighted)
# ------------------------------------------------------------
tot_df <- tot_df %>%
  mutate(
    Shop = ifelse(choice_task_outlet == "1.0", "Migros", "Coop"),
    across(matches("_quantity$"), as.numeric)
  ) %>%
  rowwise() %>%
  mutate(
    # Chicken (incl. marinated where available)
    Chicken_Cost = sum(
      PG_Chicken_quantity * (prices[[Shop]][["PG_Chicken"]] %||% 0),
      QP_Chicken_quantity * (prices[[Shop]][["QP_Chicken"]] %||% 0),
      NP_Chicken_quantity * (prices[[Shop]][["NP_Chicken"]] %||% 0),
      PG_ChickenMarinated_quantity * (prices[[Shop]][["PG_ChickenMarinated"]] %||% 0),
      Bell_ChickenMarinated_quantity * (prices[[Shop]][["Bell_ChickenMarinated"]] %||% 0),
      NP_ChickenMarinated_quantity * (prices[[Shop]][["NP_ChickenMarinated"]] %||% 0),
      na.rm = TRUE
    ),
    # Beef
    Beef_Cost = sum(
      PG_Beef_quantity * (prices[[Shop]][["PG_Beef"]] %||% 0),
      QP_Beef_quantity * (prices[[Shop]][["QP_Beef"]] %||% 0),
      NP_Beef_quantity * (prices[[Shop]][["NP_Beef"]] %||% 0),
      na.rm = TRUE
    ),
    # PBMS (Green Mountain + Planted)
    GM_Cost = sum(
      GM_Chicken_quantity * (prices[[Shop]][["GM_Chicken"]] %||% 0),
      GM_Beef_quantity * (prices[[Shop]][["GM_Beef"]] %||% 0),
      GM_ChickenMarinated_quantity * (prices[[Shop]][["GM_ChickenMarinated"]] %||% 0),
      na.rm = TRUE
    ),
    Planted_Cost = sum(
      Planted_Kebab_quantity   * (prices[[Shop]][["Planted_Kebab"]] %||% 0),
      Planted_Guggeli_quantity * (prices[[Shop]][["Planted_Guggeli"]] %||% 0),
      Planted_Skewers_quantity * (prices[[Shop]][["Planted_Skewers"]] %||% 0),
      na.rm = TRUE
    ),
    PBMS_Cost = GM_Cost + Planted_Cost,
    # Tofu
    Tofu_Cost = sum(
      Tofu_PG_quantity    * (prices[[Shop]][["Tofu_PG"]] %||% 0),
      Tofu_Karma_quantity * (prices[[Shop]][["Tofu_Karma"]] %||% 0),
      Tofu_Bio_quantity   * (prices[[Shop]][["Tofu_Bio"]] %||% 0),
      na.rm = TRUE
    ),
    # Total
    Total_Cost = Chicken_Cost + Beef_Cost + PBMS_Cost + Tofu_Cost
  ) %>%
  ungroup() %>%
  mutate(
    Chicken_Share = ifelse(Total_Cost > 0, Chicken_Cost / Total_Cost, 0),
    Beef_Share    = ifelse(Total_Cost > 0, Beef_Cost    / Total_Cost, 0),
    PBMS_Share    = ifelse(Total_Cost > 0, PBMS_Cost    / Total_Cost, 0),
    Planted_Share = ifelse(Total_Cost > 0, Planted_Cost / Total_Cost, 0),
    GM_Share      = ifelse(Total_Cost > 0, GM_Cost      / Total_Cost, 0),
    Tofu_Share    = ifelse(Total_Cost > 0, Tofu_Cost    / Total_Cost, 0),
    Share_Meat       = Chicken_Share + Beef_Share,
    Share_PlantBased = PBMS_Share,
    Share_Tofu       = Tofu_Share,
    Share_Check      = Share_Meat + Share_PlantBased + Share_Tofu
  )

# ------------------------------------------------------------
# 2) REGRESSIONS - aggregated
# ------------------------------------------------------------
dv_control_choice_meat  <- lm_robust(Share_Meat ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_meat)
dv_control_choice_pbms  <- lm_robust(Share_PlantBased ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_pbms)
dv_control_choice_tofu  <- lm_robust(Share_Tofu ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_tofu)

# Export regression table to LaTeX - Choice share variables
modelsummary(list(
  "DV4a: Meat share" = dv_control_choice_meat,
  "DV4b: PBMS share" = dv_control_choice_pbms,
  "DV4c: Tofu share" = dv_control_choice_tofu
),
  output = "latex",
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",
  gof_omit = "R2|Adj|F",
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results IIIa (stated choice shares - aggregated)",
  escape = FALSE
)

# ------------------------------------------------------------
# 2) REGRESSIONS - disaggregated
# ------------------------------------------------------------
dv_control_choice_chicken  <- lm_robust(Chicken_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_chicken)
dv_control_choice_beef  <- lm_robust(Beef_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_beef)
dv_control_choice_planted  <- lm_robust(Planted_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_planted)
dv_control_choice_gm  <- lm_robust(GM_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_gm)
dv_control_choice_tofu  <- lm_robust(Share_Tofu ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_tofu)

# Export regression table to LaTeX - Choice share variables
modelsummary(list(
  "DV4a-1: Chicken share" = dv_control_choice_chicken,
  "DV4a-2: Beef share" = dv_control_choice_beef,
  "DV4b-1: Planted share" = dv_control_choice_planted,
  "DV4b-2: Green Mountain share" = dv_control_choice_gm,
  "DV4c: Tofu share" = dv_control_choice_tofu
),
  output = "latex",
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",
  gof_omit = "R2|Adj|F",
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results IIIb (stated choice shares - disaggregated)",
  escape = FALSE
)
```

### P-adjustment
```{r}
p_values_main_analysis <- c(
  0.000000000000275,
  0.000000000000001,
  0.000000000074941,
  0.000000020305014,
  0.002334096739299,
  0.000000128054438,
  0.000861097728784,
  0.004632016824108,
  0.029153363946191,
  0.029153363946191,
  0.376746512741624,
  0.107640575826947,
  0.726132505200671
)

# Benjamini-Hochberg adjustment
p.adjust(p_values_main_analysis, method = "BH")
```


## Main results with standardized DVs
```{r}
dv_control_quality_index_std <- lm_robust(scale(quality_index) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_quality_index_std)

dv_control_health_index_std <- lm_robust(scale(health_index) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_health_index_std)

dv_control_convenience_index_std <- lm_robust(scale(convenience_index) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_convenience_index_std)

dv_control_environment_index_std <- lm_robust(scale(environment_index) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_environment_index_std)

dv_control_animal_index_std <- lm_robust(scale(animal_index) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_animal_index_std)

dv_control_price_index_std <- lm_robust(scale(price_index) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_price_index_std)


dv_control_sub_intentions_1_std <- lm_robust(scale(sub_intentions_1) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_sub_intentions_1_std)

dv_control_sub_intentions_2_std <- lm_robust(scale(sub_intentions_2) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_sub_intentions_2_std)

dv_control_meat_intentions_1_std <- lm_robust(scale(meat_intentions_1) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_meat_intentions_1_std)

dv_control_meat_intentions_2_std <- lm_robust(scale(meat_intentions_2) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_meat_intentions_2_std)

dv_control_choice_meat_std  <- lm_robust(scale(Share_Meat)  ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_meat_std)
dv_control_choice_pbms_std  <- lm_robust(scale(Share_PlantBased)  ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_pbms_std)
dv_control_choice_tofu_std  <- lm_robust(scale(Share_Tofu)  ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_tofu_std)
dv_control_choice_chicken_std  <- lm_robust(scale(Chicken_Share) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_chicken_std)
dv_control_choice_beef_std  <- lm_robust(scale(Beef_Share) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_beef_std)
dv_control_choice_planted_std  <- lm_robust(scale(Planted_Share) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_planted_std)
dv_control_choice_gm_std  <- lm_robust(scale(GM_Share) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_choice_gm_std)

dv_control_policy_support_1_std <- lm_robust(scale(policy_support_1) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_policy_support_1_std)

dv_control_policy_support_2_std <- lm_robust(scale(policy_support_2) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_policy_support_2_std)

dv_control_policy_support_3_std <- lm_robust(scale(policy_support_3) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_control_policy_support_3_std)


dv_effsize_policy_support_5_1_std <- lm(scale(policy_support_5_1) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_5_1_std)

dv_effsize_policy_support_5_2_std <- lm(scale(policy_support_5_2) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_5_2_std)

dv_effsize_policy_support_5_3_std <- lm(scale(policy_support_5_3) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_5_3_std)

dv_effsize_policy_support_5_4_std <- lm(scale(policy_support_5_4) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_5_4_std)

dv_effsize_policy_support_5_5_std <- lm(scale(policy_support_5_5) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_5_5_std)


dv_effsize_policy_support_6_1_std <- lm(scale(policy_support_6_1) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_6_1_std)

dv_effsize_policy_support_6_2_std <- lm(scale(policy_support_6_2) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_6_2_std)

dv_effsize_policy_support_6_3_std <- lm(scale(policy_support_6_3) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_6_3_std)

dv_effsize_policy_support_6_4_std <- lm(scale(policy_support_6_4) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_6_4_std)

dv_effsize_policy_support_6_5_std <- lm(scale(policy_support_6_5) ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df); summary(dv_effsize_policy_support_6_5_std)

# =========================
# Main regression results I (PBMS evaluation indexes) — standardized DVs
# =========================
modelsummary(list(
  "DV1a: Quality" = dv_control_quality_index_std,
  "DV1b: Health" = dv_control_health_index_std,
  "DV1c: Convenience" = dv_control_convenience_index_std
),
output = "latex",
stars = c('*' = .05, '**' = .01, '***' = .001),
statistic = "std.error",
gof_omit = "R2|Adj|F",
notes = list("Robust SEs in parentheses.",
             "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
title = "Main regression results I (PBMS evaluation indexes; standardized outcomes)",
escape = FALSE
)

modelsummary(list(
  "DV1d: Environment and Climate" = dv_control_environment_index_std,
  "DV1e: Animal welfare" = dv_control_animal_index_std,
  "DV1f: Price" = dv_control_price_index_std
),
output = "latex",
stars = c('*' = .05, '**' = .01, '***' = .001),
statistic = "std.error",
gof_omit = "R2|Adj|F",
notes = list("Robust SEs in parentheses.",
             "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
title = "Main regression results I (PBMS evaluation indexes; standardized outcomes — continued)",
escape = FALSE
)

# =========================
# Main regression results II (intention variables) — standardized DVs
# =========================
modelsummary(list(
  "DV2a: Reduce meat consumption in the future" = dv_control_meat_intentions_1_std,
  "DV2b: Reduce meat consumption in the next two weeks" = dv_control_meat_intentions_2_std,
  "DV3a: Purchase PBMS in the future" = dv_control_sub_intentions_1_std,
  "DV3b: Purchase PBMS in the next two weeks" = dv_control_sub_intentions_2_std
),
output = "latex",
stars = c('*' = .05, '**' = .01, '***' = .001),
statistic = "std.error",
gof_omit = "R2|Adj|F",
notes = list("Robust SEs in parentheses.",
             "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
title = "Main regression results II (intention variables; standardized outcomes)",
escape = FALSE
)

# =========================
# Main regression results III (stated choice shares) — standardized DVs
# =========================
modelsummary(list(
  "DV4a: Meat share" = dv_control_choice_meat_std,
  "DV4b: PBMS share" = dv_control_choice_pbms_std,
  "DV4c: Tofu share" = dv_control_choice_tofu_std
),
  output = "latex",
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",
  gof_omit = "R2|Adj|F",
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results III (stated choice shares; standardized outcomes)",
  escape = FALSE
)

# =========================
# Main regression results IV (general policy support) — standardized DVs
# =========================
modelsummary(list(
  "DV5a: Policies expanding plant-based food production/consumption" = dv_control_policy_support_1_std,
  "DV5b: Policies promoting PBMS" = dv_control_policy_support_2_std,
  "DV5c: Policies promoting meat reduction in CH" = dv_control_policy_support_3_std
),
output = "latex",
stars = c('*' = .05, '**' = .01, '***' = .001),
statistic = "std.error",
gof_omit = "R2|Adj|F",
notes = list("Robust SEs in parentheses.",
             "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
title = "Main regression results IV (general policy support; standardized outcomes)",
escape = FALSE
)

# =========================
# Further analyses: specific policy support (PBMS-promoting) — standardized DVs
# =========================
modelsummary(list(
  "DV6a: Lower taxes on PBMS" = dv_effsize_policy_support_5_1_std,
  "DV6b: PBMS in half of meals in public cafeterias" = dv_effsize_policy_support_5_2_std,
  "DV6c: Increased funding for research" = dv_effsize_policy_support_5_3_std
),
output = "latex",
stars = c('*' = .05, '**' = .01, '***' = .001),
statistic = "std.error",
gof_omit = "R2|Adj|F",
notes = list("Standard errors in parentheses.",
             "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
title = "Further analyses (specific policy support for measures promoting PBMS; standardized outcomes)",
escape = FALSE
)

modelsummary(list(
  "DV6d: Financial support for farmers" = dv_effsize_policy_support_5_4_std,
  "DV6e: Reduction of tariffs on PBMS raw material imports" = dv_effsize_policy_support_5_5_std
),
output = "latex",
stars = c('*' = .05, '**' = .01, '***' = .001),
statistic = "std.error",
gof_omit = "R2|Adj|F",
notes = list("Standard errors in parentheses.",
             "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
title = "Further analyses (specific policy support for measures promoting PBMS; standardized outcomes — continued)",
escape = FALSE
)

# =========================
# Further analyses: specific policy support (meat-reduction) — standardized DVs
# =========================
modelsummary(list(
  "DV7a: Higher taxes on meat products" = dv_effsize_policy_support_6_1_std,
  "DV7b: Half of meals in public cafeterias meat free" = dv_effsize_policy_support_6_2_std,
  "DV7c: Reduction of funding for research" = dv_effsize_policy_support_6_3_std
),
output = "latex",
stars = c('*' = .05, '**' = .01, '***' = .001),
statistic = "std.error",
gof_omit = "R2|Adj|F",
notes = list("Standard errors in parentheses.",
             "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
title = "Further analyses (specific policy support for measures promoting reduction of meat consumption; standardized outcomes)",
escape = FALSE
)

modelsummary(list(
  "DV7d: Less financial support for farmers producing meat" = dv_effsize_policy_support_6_4_std,
  "DV7e: Increase tariffs on animal feed imports" = dv_effsize_policy_support_6_5_std
),
output = "latex",
stars = c('*' = .05, '**' = .01, '***' = .001),
statistic = "std.error",
gof_omit = "R2|Adj|F",
notes = list("Standard errors in parentheses.",
             "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
title = "Further analyses (specific policy support for measures promoting reduction of meat consumption; standardized outcomes — continued)",
escape = FALSE
)
```


## Other Regression Models with Controls
```{r}
dv_control_policy_support_4 <- lm_robust(policy_support_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_sub_emotions <- lm_robust(sub_emotions ~ group + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_general_opinion <- lm_robust(sub_general_opinion ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_sub_purchase_reason1_1 <- lm_robust(sub_purchase_reason1_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason1_2 <- lm_robust(sub_purchase_reason1_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason1_3 <- lm_robust(sub_purchase_reason1_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason1_4 <- lm_robust(sub_purchase_reason1_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason1_5 <- lm_robust(sub_purchase_reason1_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason1_6 <- lm_robust(sub_purchase_reason1_6 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason1_7 <- lm_robust(sub_purchase_reason1_7 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)


dv_control_sub_purchase_reason2_1 <- lm_robust(sub_purchase_reason2_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason2_2 <- lm_robust(sub_purchase_reason2_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason2_3 <- lm_robust(sub_purchase_reason2_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason2_4 <- lm_robust(sub_purchase_reason2_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_sub_purchase_reason3_1 <- lm_robust(sub_purchase_reason3_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason3_2 <- lm_robust(sub_purchase_reason3_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason3_3 <- lm_robust(sub_purchase_reason3_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason3_4 <- lm_robust(sub_purchase_reason3_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_purchase_reason3_5 <- lm_robust(sub_purchase_reason3_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_sub_knowledge_1 <- lm_robust(sub_knowledge_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_sub_knowledge_2 <- lm_robust(sub_knowledge_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_sub_brand_awareness_1 <- lm_robust(sub_brand_awareness_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_social_norms_sub_1 <- lm_robust(social_norms_sub_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_social_norms_sub_2 <- lm_robust(social_norms_sub_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_social_norms_sub_3 <- lm_robust(social_norms_sub_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_social_norms_sub_4 <- lm_robust(social_norms_sub_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_social_norms_sub_5 <- lm_robust(social_norms_sub_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_control_social_norms_meat_1 <- lm_robust(social_norms_meat_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_social_norms_meat_2 <- lm_robust(social_norms_meat_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_social_norms_meat_3 <- lm_robust(social_norms_meat_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_social_norms_meat_4 <- lm_robust(social_norms_meat_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_social_norms_meat_5 <- lm_robust(social_norms_meat_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
```

# Further Analyses
## Product Subgroup Effects with Controls
### Kebab
```{r}
# Product effects with Controls
dv_sub1_control_quality_index <- lm_robust(quality_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_health_index <- lm_robust(health_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_convenience_index <- lm_robust(convenience_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_environment_index <- lm_robust(environment_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_animal_index <- lm_robust(animal_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_price_index <- lm_robust(price_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)

dv_sub1_control_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)

dv_sub1_control_policy_support_1 <- lm_robust(policy_support_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_2 <- lm_robust(policy_support_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_3 <- lm_robust(policy_support_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_4 <- lm_robust(policy_support_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)

dv_sub1_control_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)

dv_sub1_control_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)
dv_sub1_control_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = kebab_df)

# Export regression table to LaTeX - Intention Variables
modelsummary(list("DV2a: Reduce meat consumption in the future" = dv_sub1_control_meat_intentions_1,
                  "DV2b: Reduce meat consumption in the next two weeks" = dv_sub1_control_meat_intentions_2,
                  "DV3a: Purchase PBMS in the future" = dv_sub1_control_sub_intentions_1,
                  "DV3b: Purchase PBMS in the next two weeks" = dv_sub1_control_sub_intentions_2),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results II (intention variables – kebab vs. control)",
  escape = FALSE
)

# Export regression table to LaTeX - Policy Support Variables
modelsummary(list("DV5a: Policies expanding plant-based food production/consumption" = dv_sub1_control_policy_support_1,
                  "DV5b: Policies promoting PBMS" = dv_sub1_control_policy_support_2,
                  "DV5c: Policies promoting meat reduction in CH" = dv_sub1_control_policy_support_3),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results III (general policy support variables – kebab vs. control)",
  escape = FALSE
)

# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1a: Quality" = dv_sub1_control_quality_index,
                  "DV1b: Health" = dv_sub1_control_health_index,
                  "DV1c: Convenience" = dv_sub1_control_convenience_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results I (PBMS evaluation indexes – kebab vs. control)",
  escape = FALSE
)
  
  # Export regression table to LaTeX - Index Variables
  modelsummary(list("DV1d: Environment and Climate" = dv_sub1_control_environment_index,
                  "DV1e: Animal welfare" = dv_sub1_control_animal_index,
                  "DV1f: Price" =dv_sub1_control_price_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results I (PBMS evaluation indexes – kebab vs. control - continued)",
  escape = FALSE
)
```


### Guggeli
```{r}
# Product effects with Controls
dv_sub2_control_quality_index <- lm_robust(quality_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_health_index <- lm_robust(health_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_convenience_index <- lm_robust(convenience_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_environment_index <- lm_robust(environment_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_animal_index <- lm_robust(animal_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_price_index <- lm_robust(price_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)

dv_sub2_control_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)

dv_sub2_control_policy_support_1 <- lm_robust(policy_support_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_2 <- lm_robust(policy_support_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_3 <- lm_robust(policy_support_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_4 <- lm_robust(policy_support_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)

dv_sub2_control_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)

dv_sub2_control_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)
dv_sub2_control_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = guggeli_df)

# Export regression table to LaTeX - Intention Variables
modelsummary(list("DV2a: Reduce meat consumption in the future" = dv_sub2_control_meat_intentions_1,
                  "DV2b: Reduce meat consumption in the next two weeks" = dv_sub2_control_meat_intentions_2,
                  "DV3a: Purchase PBMS in the future" = dv_sub2_control_sub_intentions_1,
                  "DV3b: Purchase PBMS in the next two weeks" = dv_sub2_control_sub_intentions_2),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results II (intention variables – guggeli vs. control)",
  escape = FALSE
)

# Export regression table to LaTeX - Policy Support Variables
modelsummary(list("DV5a: Policies expanding plant-based food production/consumption" = dv_sub2_control_policy_support_1,
                  "DV5b: Policies promoting PBMS" = dv_sub2_control_policy_support_2,
                  "DV5c: Policies promoting meat reduction in CH" = dv_sub2_control_policy_support_3),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results III (general policy support variables – guggeli vs. control)",
  escape = FALSE
)

# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1a: Quality" = dv_sub2_control_quality_index,
                  "DV1b: Health" = dv_sub2_control_health_index,
                  "DV1c: Convenience" = dv_sub2_control_convenience_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results I (PBMS evaluation indexes – guggeli vs. control)",
  escape = FALSE
)
  
  # Export regression table to LaTeX - Index Variables
  modelsummary(list("DV1d: Environment and Climate" = dv_sub2_control_environment_index,
                  "DV1e: Animal welfare" = dv_sub2_control_animal_index,
                  "DV1f: Price" =dv_sub2_control_price_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results I (PBMS evaluation indexes – guggeli vs. control - continued)",
  escape = FALSE
)
```

### Skewer
```{r}
# Product effects with Controls
dv_sub3_control_quality_index <- lm_robust(quality_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_health_index <- lm_robust(health_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_convenience_index <- lm_robust(convenience_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_environment_index <- lm_robust(environment_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_animal_index <- lm_robust(animal_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_price_index <- lm_robust(price_index ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)

dv_sub3_control_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)

dv_sub3_control_policy_support_1 <- lm_robust(policy_support_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_2 <- lm_robust(policy_support_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_3 <- lm_robust(policy_support_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_4 <- lm_robust(policy_support_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)

dv_sub3_control_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)

dv_sub3_control_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)
dv_sub3_control_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = skewer_df)

# Export regression table to LaTeX - Intention Variables
modelsummary(list("DV2a: Reduce meat consumption in the future" = dv_sub3_control_meat_intentions_1,
                  "DV2b: Reduce meat consumption in the next two weeks" = dv_sub3_control_meat_intentions_2,
                  "DV3a: Purchase PBMS in the future" = dv_sub3_control_sub_intentions_1,
                  "DV3b: Purchase PBMS in the next two weeks" = dv_sub3_control_sub_intentions_2),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results II (intention variables – skewer vs. control)",
  escape = FALSE
)

# Export regression table to LaTeX - Policy Support Variables
modelsummary(list("DV5a: Policies expanding plant-based food production/consumption" = dv_sub3_control_policy_support_1,
                  "DV5b: Policies promoting PBMS" = dv_sub3_control_policy_support_2,
                  "DV5c: Policies promoting meat reduction in CH" = dv_sub3_control_policy_support_3),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results III (general policy support variables – skewer vs. control)",
  escape = FALSE
)

# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1a: Quality" = dv_sub3_control_quality_index,
                  "DV1b: Health" = dv_sub3_control_health_index,
                  "DV1c: Convenience" = dv_sub3_control_convenience_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results I (PBMS evaluation indexes – skewer vs. control)",
  escape = FALSE
)
  
  # Export regression table to LaTeX - Index Variables
  modelsummary(list("DV1d: Environment and Climate" = dv_sub3_control_environment_index,
                  "DV1e: Animal welfare" = dv_sub3_control_animal_index,
                  "DV1f: Price" =dv_sub3_control_price_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory regression results I (PBMS evaluation indexes – skewer vs. control - continued)",
  escape = FALSE
)
```

## Interaction Effects
### Group x Meat Consumption
Meat consumption variable numeric coding:
```{r}
dv_int_control_meat_quality_index <- lm_robust(quality_index ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_health_index <- lm_robust(health_index ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_convenience_index <- lm_robust(convenience_index ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_environment_index <- lm_robust(environment_index ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_animal_index <- lm_robust(animal_index ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_price_index <- lm_robust(price_index ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_meat_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_meat_policy_support_1 <- lm_robust(policy_support_1 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_2 <- lm_robust(policy_support_2 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_3 <- lm_robust(policy_support_3 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_meat_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_meat_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_meat_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group * meat_consumption + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1a: Quality" = dv_int_control_meat_quality_index,
                  "DV1b: Health" = dv_int_control_meat_health_index,
                  "DV1c: Convenience" = dv_int_control_meat_convenience_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results I (PBMS evaluation indexes)",
  escape = FALSE
)
  
# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1d: Environment and Climate" = dv_int_control_meat_environment_index,
                  "DV1e: Animal welfare" = dv_int_control_meat_animal_index,
                  "DV1f: Price" = dv_int_control_meat_price_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results I (PBMS evaluation indexes - continued)",
  escape = FALSE
)
  
# Export regression table to LaTeX - Intention Variables
modelsummary(list("DV2a: Reduce meat consumption in the future" = dv_int_control_meat_meat_intentions_1,
                  "DV2b: Reduce meat consumption in the next two weeks" = dv_int_control_meat_meat_intentions_2,
                  "DV3a: Purchase PBMS in the future" = dv_int_control_meat_sub_intentions_1,
                  "DV3a: Purchase PBMS in the next two weeks" = dv_int_control_meat_sub_intentions_2),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results II (intention variables)",
  escape = FALSE
)

# Export regression table to LaTeX - Policy Support Variables
modelsummary(list("DV5a: Policies expanding plant-based food production/consumption" = dv_int_control_meat_policy_support_1,
                  "DV5b: Policies promoting PBMS" = dv_int_control_meat_policy_support_2,
                  "DV5c: Policies promoting meat reduction in CH" = dv_int_control_meat_policy_support_3),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results III (general policy support)",
  escape = FALSE
)

```


Meat consumption as factor variable (baseline low meat):
```{r}
dv_intf_control_meat_quality_index <- lm_robust(quality_index ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_health_index <- lm_robust(health_index ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_convenience_index <- lm_robust(convenience_index ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_environment_index <- lm_robust(environment_index ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_animal_index <- lm_robust(animal_index ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_price_index <- lm_robust(price_index ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_meat_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_meat_policy_support_1 <- lm_robust(policy_support_1 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_2 <- lm_robust(policy_support_2 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_3 <- lm_robust(policy_support_3 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_meat_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_meat_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_meat_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group * meat_consumption_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption_f + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
```


### Group x Substitute Consumption
Substitute consumption variable numeric coding:
```{r}
dv_int_control_subst_quality_index <- lm_robust(quality_index ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_health_index <- lm_robust(health_index ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_convenience_index <- lm_robust(convenience_index ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_environment_index <- lm_robust(environment_index ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_animal_index <- lm_robust(animal_index ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_price_index <- lm_robust(price_index ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_subst_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_subst_policy_support_1 <- lm_robust(policy_support_1 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_2 <- lm_robust(policy_support_2 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_3 <- lm_robust(policy_support_3 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_subst_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_subst_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subst_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group * substitute_consum + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1a: Quality" = dv_int_control_subst_quality_index,
                  "DV1b: Health" = dv_int_control_subst_health_index,
                  "DV1c: Convenience" = dv_int_control_subst_convenience_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results I (PBMS evaluation indexes)",
  escape = FALSE
)
  
# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1d: Environment and Climate" = dv_int_control_subst_environment_index,
                  "DV1e: Animal welfare" = dv_int_control_subst_animal_index,
                  "DV1f: Price" = dv_int_control_subst_price_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results I (PBMS evaluation indexes - continued)",
  escape = FALSE
)
  
# Export regression table to LaTeX - Intention Variables
modelsummary(list("DV2a: Reduce meat consumption in the future" = dv_int_control_subst_meat_intentions_1,
                  "DV2b: Reduce meat consumption in the next two weeks" = dv_int_control_subst_meat_intentions_2,
                  "DV3a: Purchase PBMS in the future" = dv_int_control_subst_sub_intentions_1,
                  "DV3a: Purchase PBMS in the next two weeks" = dv_int_control_subst_sub_intentions_2),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results II (intention variables)",
  escape = FALSE
)

# Export regression table to LaTeX - Policy Support Variables
modelsummary(list("DV5a: Policies expanding plant-based food production/consumption" = dv_int_control_subst_policy_support_1,
                  "DV5b: Policies promoting PBMS" = dv_int_control_subst_policy_support_2,
                  "DV5c: Policies promoting meat reduction in CH" = dv_int_control_subst_policy_support_3),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results III (general policy support)",
  escape = FALSE
)

```

Substitute consumption as factor variable (baseline no prior substitute consumption):
```{r}
dv_intf_control_subst_quality_index <- lm_robust(quality_index ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_health_index <- lm_robust(health_index ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_convenience_index <- lm_robust(convenience_index ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_environment_index <- lm_robust(environment_index ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_animal_index <- lm_robust(animal_index ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_price_index <- lm_robust(price_index ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_subst_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_subst_policy_support_1 <- lm_robust(policy_support_1 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_2 <- lm_robust(policy_support_2 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_3 <- lm_robust(policy_support_3 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_subst_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_subst_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subst_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group * substitute_consum_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum_f + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
```


### Group x Substitute Opinion
Substitute general opinion variable numeric coding (baseline prior negative substitute opinion):
```{r}
dv_int_control_subopinion_quality_index <- lm_robust(quality_index ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_health_index <- lm_robust(health_index ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_convenience_index <- lm_robust(convenience_index ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_environment_index <- lm_robust(environment_index ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_animal_index <- lm_robust(animal_index ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_price_index <- lm_robust(price_index ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_subopinion_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_subopinion_policy_support_1 <- lm_robust(policy_support_1 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_2 <- lm_robust(policy_support_2 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_3 <- lm_robust(policy_support_3 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_subopinion_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group * sub_general_opinion_w1+ age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

dv_int_control_subopinion_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_int_control_subopinion_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group * sub_general_opinion_w1 + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1a: Quality" = dv_int_control_subopinion_quality_index,
                  "DV1b: Health" = dv_int_control_subopinion_health_index,
                  "DV1c: Convenience" = dv_int_control_subopinion_convenience_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results I (PBMS evaluation indexes)",
  escape = FALSE
)
  
# Export regression table to LaTeX - Index variables
  modelsummary(list("DV1d: Environment and Climate" = dv_int_control_subopinion_environment_index,
                  "DV1e: Animal welfare" = dv_int_control_subopinion_animal_index,
                  "DV1f: Price" = dv_int_control_subopinion_price_index),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results I (PBMS evaluation indexes - continued)",
  escape = FALSE
)
  
# Export regression table to LaTeX - Intention Variables
modelsummary(list("DV2a: Reduce meat consumption in the future" = dv_int_control_subopinion_meat_intentions_1,
                  "DV2b: Reduce meat consumption in the next two weeks" = dv_int_control_subopinion_meat_intentions_2,
                  "DV3a: Purchase PBMS in the future" = dv_int_control_subopinion_sub_intentions_1,
                  "DV3a: Purchase PBMS in the next two weeks" = dv_int_control_subopinion_sub_intentions_2),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results II (intention variables)",
  escape = FALSE
)

# Export regression table to LaTeX - Policy Support Variables
modelsummary(list("DV5a: Policies expanding plant-based food production/consumption" = dv_int_control_subopinion_policy_support_1,
                  "DV5b: Policies promoting PBMS" = dv_int_control_subopinion_policy_support_2,
                  "DV5c: Policies promoting meat reduction in CH" = dv_int_control_subopinion_policy_support_3),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Interaction regression results III (general policy support)",
  escape = FALSE
)
```

Substitute general opinion as factor variable:
```{r}
dv_intf_control_subopinion_quality_index <- lm_robust(quality_index ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_health_index <- lm_robust(health_index ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_convenience_index <- lm_robust(convenience_index ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_environment_index <- lm_robust(environment_index ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_animal_index <- lm_robust(animal_index ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_price_index <- lm_robust(price_index ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_subopinion_sub_intentions_1 <- lm_robust(sub_intentions_1 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_sub_intentions_2 <- lm_robust(sub_intentions_2 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_meat_intentions_1 <- lm_robust(meat_intentions_1 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_meat_intentions_2 <- lm_robust(meat_intentions_2 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_subopinion_policy_support_1 <- lm_robust(policy_support_1 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_2 <- lm_robust(policy_support_2 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_3 <- lm_robust(policy_support_3 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_subopinion_policy_support_5_1 <- lm_robust(policy_support_5_1 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_5_2 <- lm_robust(policy_support_5_2 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_5_3 <- lm_robust(policy_support_5_3 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_5_4 <- lm_robust(policy_support_5_4 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_5_5 <- lm_robust(policy_support_5_5 ~ group * sub_general_opinion_w1_f+ age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)

dv_intf_control_subopinion_policy_support_6_1 <- lm_robust(policy_support_6_1 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_6_2 <- lm_robust(policy_support_6_2 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_6_3 <- lm_robust(policy_support_6_3 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_6_4 <- lm_robust(policy_support_6_4 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
dv_intf_control_subopinion_policy_support_6_5 <- lm_robust(policy_support_6_5 ~ group * sub_general_opinion_w1_f + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1_f + sub_brand_awareness_1_w1, data = tot_df)
```

## Exploratory Stated Choice Preferences (conditional on purchase in the online shopping task)
```{r}
dv_control_stated_meat <- lm_robust(Share_Meat ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_stated_PBMS <- lm_robust(Share_PlantBased ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
dv_control_stated_tofu <- lm_robust(Tofu_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

# conditional on buying - aggregated
dv_control_stated_meat_positive <- lm_robust(Share_Meat ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1 + Total_Cost, data = subset(tot_df, tot_df$Share_Check==1))
dv_control_stated_PBMS_positive <- lm_robust(Share_PlantBased ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1 + Total_Cost, data = subset(tot_df, tot_df$Share_Check==1))
dv_control_stated_tofu_positive <- lm_robust(Tofu_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1+ Total_Cost, data = subset(tot_df, tot_df$Share_Check==1))

# conditional on buying - disaggregated
dv_control_stated_chicken_positive <- lm_robust(Chicken_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1 + Total_Cost, data = subset(tot_df, tot_df$Share_Check==1))
dv_control_stated_beef_positive <- lm_robust(Beef_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1 + Total_Cost, data = subset(tot_df, tot_df$Share_Check==1))
dv_control_stated_Planted_positive <- lm_robust(Planted_Share  ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1 + Total_Cost, data = subset(tot_df, tot_df$Share_Check==1))
dv_control_stated_GM_positive <- lm_robust(GM_Share  ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1 + Total_Cost, data = subset(tot_df, tot_df$Share_Check==1))
dv_control_stated_tofu_positive <- lm_robust(GM_Share ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1 + Total_Cost, data = subset(tot_df, tot_df$Share_Check==1))

# Export regression table to LaTeX - Stated choice aggregated
modelsummary(list("DV4a: Stated Choice Share of Meat" = dv_control_stated_meat,
                  "DV4b: Stated Choice Share of PBMS" = dv_control_stated_PBMS,
                  "DV4c: Stated Choice Share of Tofu" = dv_control_stated_tofu),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results II (intention variables)",
  escape = FALSE
)

# Export regression table to LaTeX - Stated choice aggregated positive only
modelsummary(list("DV4a: Stated Choice Share of Meat" = dv_control_stated_meat_positive,
                  "DV4b: Stated Choice Share of PBMS" = dv_control_stated_PBMS_positive,
                  "DV4c: Stated Choice Share of Tofu" = dv_control_stated_tofu_positive),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results II (intention variables)",
  escape = FALSE
)

# Export regression table to LaTeX - Stated Choice disagreggated
modelsummary(list("DV4a.1: Stated Choice Share of Chicken" = dv_control_stated_chicken_positive,
                  "DV4a.2: Stated Choice Share of Beef" = dv_control_stated_beef_positive,
                  "DV4b.1: Stated Choice Share of Planted" = dv_control_stated_Planted_positive,
                  "DV4b.2: Stated Choice Share of GM PBMS" = dv_control_stated_GM_positive,
                  "DV4c: Stated Choice Share of Tofu" = dv_control_stated_tofu),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Exploratory analysis stated choice III (intention variables)",
  escape = FALSE
)

modelsummary(list("DV4a.1: Stated Choice Share of Chicken" = dv_control_stated_chicken_positive,
                  "DV4a.2: Stated Choice Share of Beef" = dv_control_stated_beef_positive,
                  "DV4b.1: Stated Choice Share of Planted" = dv_control_stated_Planted_positive,
                  "DV4b.2: Stated Choice Share of GreenMountain" = dv_control_stated_GM_positive,
                  "DV4c: Stated Choice Share of Tofu" = dv_control_stated_tofu),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Main regression results II (intention variables)",
  escape = FALSE
)
```



# Manipulation and Robustness Checks
## Manipulation Check
```{r}
tot_df$planted_evaluation <- as.numeric(tot_df$planted_evaluation)
summary(tot_df$planted_evaluation)
tot_df$planted_frequency <- as.numeric(tot_df$planted_frequency)
summary(tot_df$planted_frequency)
mc_control_brandawareness <- lm_robust(sub_brand_awareness_1 ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
mc_control_plantedeval <- lm_robust(planted_evaluation ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)
mc_control_plantedfreq <- lm_robust(planted_frequency ~ group + age + gender + education + diet + left_right + urban_rural_true + meat_consumption + sub_general_opinion_w1 + sub_brand_awareness_1_w1, data = tot_df)

# Export regression table to LaTeX - Manipulation Check
  modelsummary(list("Planted brand awareness" = mc_control_brandawareness,
                  "Planted consumption frequency" = mc_control_plantedfreq,
                  "Planted evaluation" = mc_control_plantedeval),
             output = "latex",             # direct LaTeX output
  stars = c('*' = .05, '**' = .01, '***' = .001),
  statistic = "std.error",      # shows standard errors in ()
  gof_omit = "R2|Adj|F",        # hides unnecessary fit stats
  notes = list("Robust SEs in parentheses.",
               "+ p<0.1, * p<0.05, ** p<0.01, *** p<0.001"),
  title = "Manipulation check results",
  escape = FALSE
)
```

## Robustness Check
### Inverse probability weighting
```{r}
## Guggeli
# Subset to control + any guggeli groups
df_guggeli <- tot_df %>%
  filter(product_treated %in% c("control",
                                "guggeli_only",
                                "guggeliANDkebab",
                                "guggeliANDskewer",
                                "all_products")) %>%
  mutate(treat = ifelse(product_treated == "control", 0, 1))

# Define formula for covariates
form <- treat ~ age + gender + education + diet + left_right + urban_rural_true +
                meat_consumption + substitute_consum + sub_general_opinion_w1 +
                sub_brand_awareness_1_w1

# Estimate weights
w_guggeli <- weightit(form, data = df_guggeli, method = "ps", estimand = "ATE")

# Balance plot
love.plot(bal.tab(w_guggeli), thresholds = 0.1, abs = TRUE, var.order = "unadjusted", title = "Balance Plot: Guggeli vs Control")

# run IPW and create plot
love.plot(
  treat ~ age + gender + education + diet + left_right + urban_rural_true +
           meat_consumption + substitute_consum + sub_general_opinion_w1 + sub_brand_awareness_1_w1,
  data    = df_guggeli,
  weights = w_guggeli$weights,
  thresholds = 0.1, abs = TRUE, var.order = "unadjusted", stars = "raw"
)

## Kebab
# Subset to control + any kebab groups
df_kebab <- tot_df %>%
  filter(product_treated %in% c("control",
                                "kebab_only",
                                "guggeliANDkebab",
                                "kebabANDskewer",
                                "all_products")) %>%
  mutate(treat = ifelse(product_treated == "control", 0, 1))

# Define covariate formula
form <- treat ~ age + gender + education + diet + left_right + urban_rural_true +
                meat_consumption + substitute_consum + sub_general_opinion_w1 +
                sub_brand_awareness_1_w1

# Estimate weights
w_kebab <- weightit(form, data = df_kebab, method = "ps", estimand = "ATE")

# Balance plot
love.plot(bal.tab(w_kebab), thresholds = 0.1, abs = TRUE, var.order = "unadjusted", title = "Balance Plot: Kebab vs Control")

## Skewers

# Subset to control + any skewer groups
df_skewer <- tot_df %>%
  filter(product_treated %in% c("control",
                                "skewer_only",
                                "guggeliANDskewer",
                                "kebabANDskewer",
                                "all_products")) %>%
  mutate(treat = ifelse(product_treated == "control", 0, 1))

# Define covariate formula
form <- treat ~ age + gender + education + diet + left_right + urban_rural_true +
                meat_consumption + substitute_consum + sub_general_opinion_w1 +
                sub_brand_awareness_1_w1

# Estimate weights
w_skewer <- weightit(form, data = df_skewer, method = "ps", estimand = "ATE")

# Balance plot
love.plot(bal.tab(w_skewer), thresholds = 0.1, abs = TRUE, var.order = "unadjusted", title = "Balance Plot: Skewer vs Control")

# Display the combined plot
plot_guggeli <- love.plot(bal.tab(w_guggeli), thresholds=0.1, abs=TRUE,
                          var.order="unadjusted", title="Guggeli vs Control")
plot_kebab   <- love.plot(bal.tab(w_kebab), thresholds=0.1, abs=TRUE,
                          var.order="unadjusted", title="Kebab vs Control")
plot_skewer  <- love.plot(bal.tab(w_skewer), thresholds=0.1, abs=TRUE,
                          var.order="unadjusted", title="Skewer vs Control")

combined_plot <- plot_guggeli / plot_kebab / plot_skewer  # stack vertically
combined_plot





```

### Models and plots with IPW
```{r}
# Covariates used in the propensity score and in the outcome model
covs <- c(
  "age","gender","education","diet","left_right","urban_rural_true",
  "meat_consumption","substitute_consum","sub_general_opinion_w1",
  "sub_brand_awareness_1_w1"
)

# Which arms include each product (adjust if your labels differ)
product_sets <- list(
  guggeli = c("guggeli_only","guggeliANDkebab","guggeliANDskewer","all_products"),
  kebab   = c("kebab_only","guggeliANDkebab","kebabANDskewer","all_products"),
  skewer  = c("skewer_only","guggeliANDskewer","kebabANDskewer","all_products")
)

# Desired y-axis order (top to bottom) — ggplot puts the *last* level on top,
# so we pass the reversed order to scale_y_discrete below.
y_limits <- c("Skewer","Kebab","Guggeli","Overall Effect")

# =========================
# HELPERS
# =========================
# Build subset for each contrast and define treat (1 = any treated, 0 = control)
make_contrast_data <- function(df, contrast = c("overall","guggeli","kebab","skewer")) {
  contrast <- match.arg(contrast)
  df <- df %>% mutate(treat = as.integer(product_treated != "control"))

  if (contrast == "overall") {
    keep <- c("control","guggeli_only","kebab_only","skewer_only",
              "guggeliANDkebab","guggeliANDskewer","kebabANDskewer","all_products")
  } else {
    keep <- c("control", product_sets[[contrast]])
  }
  df %>% filter(product_treated %in% keep)
}

# Fit PS weights (ATE) + weighted outcome model (doubly-robust if add_controls = TRUE)
fit_weighted_te <- function(data, dv, covs, add_controls = TRUE) {
  f_ps <- as.formula(paste("treat ~", paste(covs, collapse = " + ")))
  w <- WeightIt::weightit(f_ps, data = data, method = "ps", estimand = "ATE")

  rhs <- if (add_controls) paste("treat +", paste(covs, collapse = " + ")) else "treat"
  f_y <- as.formula(paste(dv, "~", rhs))
  fit <- estimatr::lm_robust(f_y, data = data, weights = w$weights, se_type = "HC2")

  list(weights = w, fit = fit)
}

# Run IPW for all contrasts (Overall/Guggeli/Kebab/Skewer) × DVs, return tidy table
run_all_weighted <- function(df, dv_named_vec, covs, add_controls = TRUE) {
  contrasts <- c("overall","guggeli","kebab","skewer")
  out <- list()

  for (ct in contrasts) {
    d_ct <- make_contrast_data(df, ct)
    for (dv_lab in names(dv_named_vec)) {
      dv_var <- dv_named_vec[[dv_lab]]
      res <- fit_weighted_te(d_ct, dv_var, covs, add_controls)

      # Robust effective sample size extraction (optional meta)
      ess <- tryCatch(summary(res$weights)$eff.n, error = function(e) NULL)
      ess_ctrl <- NA_real_; ess_trt <- NA_real_
      if (!is.null(ess)) {
        nms <- tolower(names(ess))
        if ("control" %in% nms) ess_ctrl <- as.numeric(ess[which(nms=="control")[1]])
        if ("treated" %in% nms) ess_trt  <- as.numeric(ess[which(nms=="treated")[1]])
        if (is.na(ess_ctrl)) ess_ctrl <- as.numeric(ess[1])
        if (is.na(ess_trt))  ess_trt  <- as.numeric(ess[length(ess)])
      }

      est <- broom::tidy(res$fit, conf.int = TRUE) |>
        dplyr::filter(term == "treat") |>
        dplyr::mutate(
          contrast = ct,
          dv_label = dv_lab,
          dv       = dv_var,
          ess_ctrl = round(ess_ctrl, 2),
          ess_trt  = round(ess_trt,  2)
        )
      out[[length(out)+1]] <- est
    }
  }
  dplyr::bind_rows(out)
}

# Plotter used for *all* figures (style matches your p_reasons)
# IMPORTANT: this function reads facet order from global `dv_list`.
plot_reasons_panel <- function(tidy_tbl, title) {
  tidy_tbl$contrast <- factor(
    tidy_tbl$contrast,
    levels = c("overall","guggeli","kebab","skewer"),
    labels = c("Overall Effect","Guggeli","Kebab","Skewer")
  )
  tidy_tbl$dv_label <- factor(tidy_tbl$dv_label, levels = names(dv_list))

  ggplot(tidy_tbl, aes(x = estimate, y = contrast)) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_point(size = 2.4, color = "blue") +
    geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
                   height = 0.15, color = "blue") +
    facet_wrap(~ dv_label, ncol = 3, scales = "free_x") +
    labs(x = "Estimate", y = NULL, title = title) +
    theme_minimal(base_size = 11) +
    theme(
      legend.position = "none",
      panel.background = element_rect(fill = "grey90"),
      panel.grid.major = element_line(color = "white"),
      strip.text = element_text(face = "bold"),
      plot.title = element_text(hjust = 0.5, face = "bold")
    ) +
    scale_y_discrete(limits = y_limits)  # top-to-bottom: Overall, Guggeli, Kebab, Skewer
}

# Helper to reuse the same plotter with different facet orders (no edits to the plotter)
with_dv_levels <- function(dv_levels, expr) {
  old <- dv_list
  on.exit(dv_list <<- old, add = TRUE)
  dv_list <<- dv_levels
  force(expr)
}

# =========================
# FIGURE 1: Purchase Reasons (p_reasons)
# =========================
dv_reasons <- c(
  "Quality"        = "sub_purchase_reason1_1",
  "Health"         = "sub_purchase_reason1_2",
  "Price"          = "sub_purchase_reason1_3",
  "Convenience"    = "sub_purchase_reason1_4",
  "Environment"    = "sub_purchase_reason1_5",
  "Animal Welfare" = "sub_purchase_reason1_7"
)

tidy_reasons <- run_all_weighted(tot_df, dv_reasons, covs, add_controls = TRUE)

dv_list <- NULL 

p_reasons <- with_dv_levels(
  dv_levels = dv_reasons,
  expr = plot_reasons_panel(
    tidy_reasons,
    title = "Reasons to Purchase PBMA (Weighted IPW)"
  )
)
print(p_reasons)
# ggsave("plot_reasons_weighted.png", p_reasons, width=20, height=14, units="cm", dpi=300)

# =========================
# FIGURE 2: Intentions (p_intentions)
# =========================
dv_intentions <- c(
  "PBMA in the Future"                            = "sub_intentions_1",
  "PBMA in the next two Weeks"                    = "sub_intentions_2",
  "Reduce Meat Consumption in the Future"         = "meat_intentions_1",
  "Reduce Meat consumption in the next two Weeks" = "meat_intentions_2"
)

tidy_intentions <- run_all_weighted(tot_df, dv_intentions, covs, add_controls = TRUE)

p_intentions <- with_dv_levels(
  dv_levels = dv_intentions,
  expr = plot_reasons_panel(
    tidy_intentions,
    title = "Intentions to Consume PBMA and Meat in the Future (Weighted IPW)"
  )
)
print(p_intentions)
# ggsave("plot_intentions_weighted.png", p_intentions, width=20, height=14, units="cm", dpi=300)

# =========================
# FIGURE 3: Policy Support (p_policy)
# =========================
# Use any set you like; here’s a simple two-facet example:
dv_policy <- c(
  "Increase production and consumption of plant-based foods in Switzerland" = "policy_support_1",
  "Promotion of PBMA in Switzerland"                                        = "policy_support_2",
  "Reduction of Meat Consumption in Switzerland"                            = "policy_support_3"
  # Add more policies here if needed
)

tidy_policy <- run_all_weighted(tot_df, dv_policy, covs, add_controls = TRUE)

p_policy <- with_dv_levels(
  dv_levels = dv_policy,
  expr = plot_reasons_panel(
    tidy_policy,
    title = "Policy Support (Weighted IPW)"
  )
)
print(p_policy)

## Share
# ------------------------------------------------------------
# PRICE LISTS (per package, as in your earlier script)
# ------------------------------------------------------------
prices <- list(
  Migros = list(
    PG_Chicken = 9.40, QP_Chicken = 10.55, NP_Chicken = 11.45,
    PG_Beef = 11.20, QP_Beef = 16.00, NP_Beef = 13.50,
    Planted_Kebab = 5.95, Planted_Guggeli = 5.95, Planted_Skewers = 6.95,
    Tofu_PG = 1.55, Tofu_Karma = 3.30, Tofu_Bio = 3.50
  ),
  Coop = list(
    PG_Chicken = 7.05, QP_Chicken = 12.00, NP_Chicken = 12.30,
    PG_Beef = 11.70, QP_Beef = 15.80, NP_Beef = 13.10,
    PG_ChickenMarinated = 2.95, Bell_ChickenMarinated = 7.90, NP_ChickenMarinated = 11.10,
    GM_Chicken = 5.95, GM_Beef = 5.95, GM_ChickenMarinated = 5.95,
    Planted_Kebab = 5.95, Planted_Guggeli = 5.95, Planted_Skewers = 6.95,
    Tofu_PG = 1.60, Tofu_Karma = 2.95, Tofu_Bio = 3.95
  )
)

# ------------------------------------------------------------
# 1) CALCULATE STATED-CHOICE SHARES WITH tot_df
# ------------------------------------------------------------

tot_df <- tot_df %>%
  mutate(
    Shop = ifelse(choice_task_outlet == "1.0", "Migros", "Coop"),
    across(matches("_quantity$"), as.numeric)
  )

tot_df <- tot_df %>%
  rowwise() %>%
  mutate(
    Chicken_Cost = sum(
      PG_Chicken_quantity * prices[[Shop]][["PG_Chicken"]],
      QP_Chicken_quantity * prices[[Shop]][["QP_Chicken"]],
      NP_Chicken_quantity * prices[[Shop]][["NP_Chicken"]],
      PG_ChickenMarinated_quantity * (prices[[Shop]][["PG_ChickenMarinated"]] %||% 0),
      Bell_ChickenMarinated_quantity * (prices[[Shop]][["Bell_ChickenMarinated"]] %||% 0),
      NP_ChickenMarinated_quantity * (prices[[Shop]][["NP_ChickenMarinated"]] %||% 0),
      na.rm = TRUE
    ),
    Beef_Cost = sum(
      PG_Beef_quantity * prices[[Shop]][["PG_Beef"]],
      QP_Beef_quantity * prices[[Shop]][["QP_Beef"]],
      NP_Beef_quantity * prices[[Shop]][["NP_Beef"]],
      na.rm = TRUE
    ),
    GM_Cost = sum(
      GM_Chicken_quantity * (prices[[Shop]][["GM_Chicken"]] %||% 0),
      GM_Beef_quantity * (prices[[Shop]][["GM_Beef"]] %||% 0),
      GM_ChickenMarinated_quantity * (prices[[Shop]][["GM_ChickenMarinated"]] %||% 0),
      na.rm = TRUE
    ),
    Planted_Cost = sum(
      Planted_Kebab_quantity   * prices[[Shop]][["Planted_Kebab"]],
      Planted_Guggeli_quantity * prices[[Shop]][["Planted_Guggeli"]],
      Planted_Skewers_quantity * prices[[Shop]][["Planted_Skewers"]],
      na.rm = TRUE
    ),
    Tofu_Cost = sum(
      Tofu_PG_quantity    * prices[[Shop]][["Tofu_PG"]],
      Tofu_Karma_quantity * prices[[Shop]][["Tofu_Karma"]],
      Tofu_Bio_quantity   * prices[[Shop]][["Tofu_Bio"]],
      na.rm = TRUE
    ),
    Total_Cost = Chicken_Cost + Beef_Cost + GM_Cost + Planted_Cost + Tofu_Cost
  ) %>%
  ungroup() %>%
  mutate(
    Chicken_Share = ifelse(Total_Cost > 0, Chicken_Cost / Total_Cost, 0),
    Beef_Share    = ifelse(Total_Cost > 0, Beef_Cost    / Total_Cost, 0),
    GM_Share      = ifelse(Total_Cost > 0, GM_Cost      / Total_Cost, 0),
    Planted_Share = ifelse(Total_Cost > 0, Planted_Cost / Total_Cost, 0),
    Tofu_Share    = ifelse(Total_Cost > 0, Tofu_Cost    / Total_Cost, 0),
    Share         = round(Chicken_Share + Beef_Share + GM_Share + Planted_Share + Tofu_Share, 3),
    Share_Meat        = Chicken_Share + Beef_Share,
    Share_PlantBased  = GM_Share + Planted_Share
  )

# ------------------------------------------------------------
# 2) RUN IPW MODELS ON SHARES (with your existing helpers)
# ------------------------------------------------------------

dv_shares <- c(
  "Meat Share"           = "Share_Meat",
  "Planted Share"        = "Planted_Share",
  "Green Mountain Share" = "GM_Share",
  "Tofu Share"           = "Tofu_Share"
)

tidy_shares <- run_all_weighted(tot_df, dv_shares, covs, add_controls = TRUE)

p_shares <- with_dv_levels(
  dv_levels = dv_shares,
  expr = plot_reasons_panel(
    tidy_shares,
    title = "Stated-Choice: Product Shares (Weighted IPW)"
  )
)

print(p_shares)


# Save Plots
ggsave("plot_reasons_ipw.png", p_reasons,
       height = 14, width = 20, units = "cm", dpi = 300)

# Save Policy plot
ggsave("plot_policy_ipw.png", p_policy,
       height = 14, width = 20, units = "cm", dpi = 300)

# Save Intentions plot
ggsave("plot_intentions_ipw.png", p_intentions,
       height = 15, width = 15, units = "cm", dpi = 300)

# Save Stated Choices plot
ggsave("stated_choice_shares_ipw.png", p_shares,
       width = 20, height = 14, units = "cm", dpi = 300)

# Save IPW check
ggsave("ipw_check.png", combined_plot,
       height = 15, width = 15, units = "cm", dpi = 300)

```

# Plots
## Average Treatment Effect Plots
```{r}
#Summary plot for main effects: Indeces
main_effects_1 <- plot_summs(
  dv_control_quality_index, dv_control_health_index, dv_control_price_index,
  dv_control_convenience_index, dv_control_environment_index, dv_control_animal_index,
  scale = TRUE, inner_ci_level = 0.9,
  coefs = c("Experience Treatment" = "group"),
  model.names = c("Quality", "Health", "Price", "Convenience",
                  "Climate and Environment", "Animal Welfare"),
  colors = c("#1f78b4", "#cab2d6", "#999999", "#6a3d9a", "#b15928", "#ff7f00")
) +
  xlim(-0.2, 0.6) +
  scale_y_discrete(expand = expansion(mult = c(0, 0.01))) +  # <— tighten top/bottom+
  theme(
    legend.position = "bottom",
    axis.title = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.text  = element_text(size = 9),
    axis.text    = element_text(size = 10),
    axis.text.y  = element_text(size = 10),
    plot.margin  = unit(c(0, 0.25, 0, 0.25), "cm")            # <— trim outer margins
  )

ggsave("plot_indeces_ate.pdf", main_effects_1, height = 14, width = 20,
       units = "cm", dpi = 300, limitsize = FALSE)

#Summary plot for main effects: Intentions
main_effects_2 <- plot_summs(
  dv_control_sub_intentions_1, dv_control_sub_intentions_2, dv_control_meat_intentions_1,
  dv_control_meat_intentions_2,
  scale = TRUE, inner_ci_level = 0.9,
  coefs = c("Experience Treatment" = "group"),
  model.names = c("Increase PBMS \nnext 6 months", "PBMS \nnext 2 weeks", "Reduce meat \nnext 6 months", "Reduce meat \nnext 2 weeks"),
  colors = c("#b2df8a","#33a02c","#fb9a99","#e31a1c")
) +
  xlim(-0.2, 0.6)+
  theme(legend.position="bottom",axis.title=element_text(size=10), legend.title = element_text(size = 10),legend.text = element_text(size = 9), axis.text=element_text(size=10), axis.text.y = element_text(size = 10)) 

ggsave("plot_intentions_ate.pdf", main_effects_2, height = 14, width = 20, units = "cm", dpi = 300, limitsize = FALSE)

#Summary plot for main effects: Stated product choices
main_effects_choice <- plot_summs(
  dv_control_choice_meat, dv_control_choice_pbms, dv_control_choice_tofu,
  scale = TRUE, inner_ci_level = 0.9,
  coefs = c("Experience Treatment" = "group"),
  model.names = c("Meat share", "PBMS share", "Tofu share"),
  colors = c("#e31a1c", "#33a02c", "#1f78b4")
) +
  xlim(-0.2, 0.6) +
  scale_y_discrete(expand = expansion(mult = c(0, 0.01))) +
  xlab("Estimate (standardized)") +
  theme(
    legend.position = "bottom",
    axis.title = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.text  = element_text(size = 9),
    axis.text    = element_text(size = 10),
    axis.text.y  = element_text(size = 10),
    plot.margin  = unit(c(0, 0.25, 0, 0.25), "cm")
  )

ggsave("plot_choice_shares_ate.pdf", main_effects_choice,
       height = 14, width = 20, units = "cm", dpi = 300, limitsize = FALSE)

#Summary plot for main effects: Policy support PBMS
main_effects_3 <- plot_summs(
  dv_control_policy_support_5_1, dv_control_policy_support_5_2, dv_control_policy_support_5_3,
  dv_control_policy_support_5_4, dv_control_policy_support_5_5,
  scale = TRUE, inner_ci_level = 0.9,
  coefs = c("Experience Treatment" = "group"),
  model.names = c(
    "Lower taxes \non PBMS",
    "PBMS in 50% of public \ncafeteria meals",
    "Increased funding \nfor PBMS R&D",
    "Increase support for farmers \nproducing PBMS raw materials",
    "Reduction of tariffs on \nPBMS raw material imports"
  ),
  colors = c("#009E73", "#56B4E9", "#CC79A7", "#a6cee3", "#1f78b4")
) +
  xlim(-0.2, 0.6) +
  scale_y_discrete(expand = expansion(mult = c(0.03, 0.12))) +  # more space esp. on top
  coord_cartesian(clip = "off") +                               # prevent clipping
  theme(
    legend.position = "bottom",
    axis.title = element_text(size = 10),
    legend.title = element_text(size = 10),
    legend.text  = element_text(size = 8),
    legend.key.width = unit(0.6, "cm"),
    legend.spacing.x = unit(0.2, "cm"),
    axis.text    = element_text(size = 10),
    axis.text.y  = element_text(size = 10),
    plot.margin  = margin(t = 6, r = 8, b = 6, l = 8, unit = "pt")  # a bit more top/right
  )  +
  guides(color = guide_legend(nrow = 2, byrow = TRUE))   # <- key line

ggsave("plot_policy_support_1_ate.pdf", main_effects_3, height = 14, width = 20, units = "cm", dpi = 300, limitsize = FALSE)

#Summary plot for main effects: Policy support Meat
main_effects_4 <- plot_summs(
  dv_control_policy_support_3, dv_control_policy_support_6_1, dv_control_policy_support_6_2, dv_control_policy_support_6_3, dv_control_policy_support_6_4, dv_control_policy_support_6_5,
  scale = TRUE, inner_ci_level = 0.9,
  coefs = c("Experience Treatment" = "group"),
  model.names = c("General reduction of \nmeat consumption in CH", "Higher taxes on \nmeat products", "No meat in 50% of \npublic cafeteria meals", "Reduce funding for \nmeat product R&D", "Less financial support \nfor farmers producing \nmeat & animal feed", "Increase tariffs on \nanimal feed imports"),
  colors = c("darkgrey", "#E69F00", "#F0E442", "#D55E00", "darkred", "red")
) +
  xlim(-0.4, 0.4)+
  theme(legend.position="bottom",axis.title=element_text(size=10), legend.title = element_text(size = 10),legend.text = element_text(size = 9), axis.text=element_text(size=10), axis.text.y = element_text(size = 10)) 

ggsave("plot_policy_support_2_ate.pdf", main_effects_4, height = 14, width = 20, units = "cm", dpi = 300, limitsize = FALSE)

```

## Interaction Effect Plots
```{r}
# Refit model with lm() for interaction plotting
dv_int_control_subst_sub_intentions_1_lm <- lm(
  sub_intentions_1 ~ group * substitute_consum + age + gender + education + diet +
    left_right + urban_rural_true + meat_consumption + substitute_consum +
    sub_general_opinion_w1 + sub_brand_awareness_1_w1,
  data = tot_df)

#Interaction effect (Group x prior substitute experience): Substitute consumption intentions
plot_interact_intentions_subst <- interact_plot(
  dv_int_control_subst_sub_intentions_1_lm,
  pred = substitute_consum,
  modx = group,
  interval = TRUE,
  int.type = "confidence",
  int.width = 0.95,
  modx.labels = c("Control", "Experience Treatment"),
  legend.main = "Group",
  colors = c("#1f78b4", "#33a02c") 
) +
  xlab("Prior PBMS consumption") +
  ylab("Intention to increase PBMS consumption (next 6 months)") +
  theme_classic(base_size = 10) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.text = element_text(size = 9),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    axis.line = element_line(linewidth = 0.8),
    axis.ticks = element_line(linewidth = 0.8),
    panel.grid = element_blank(),       # removes all grid lines
    plot.title = element_text(size = 11, face = "bold", hjust = 0.5),
    plot.margin = margin(5, 5, 5, 5)
  )

# Save as high-quality PDF for publication
ggsave(
  "plot_interaction_subst_intentions.pdf",
  plot_interact_intentions_subst,
  height = 14, width = 20, units = "cm",
  dpi = 600, limitsize = FALSE
)

# Refit model with lm() for interaction plotting
dv_int_control_subst_meat_intentions_1_lm <- lm(
  meat_intentions_1 ~ group * substitute_consum + age + gender + education + diet +
    left_right + urban_rural_true + meat_consumption + substitute_consum +
    sub_general_opinion_w1 + sub_brand_awareness_1_w1,
  data = tot_df
)

# Interaction effect (Group x prior substitute experience): Meat reduction intentions
plot_interact_intentions_meat <- interact_plot(
  dv_int_control_subst_meat_intentions_1_lm,
  pred = substitute_consum,
  modx = group,
  interval = TRUE,
  int.type = "confidence",
  int.width = 0.95,
  modx.labels = c("Control", "Experience Treatment"),
  legend.main = "Group",
  colors = c("#1f78b4", "#e31a1c") 
) +
  xlab("Prior PBMS consumption") +
  ylab("Intention to reduce meat consumption (next 6 months)") +
  theme_classic(base_size = 10) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.text = element_text(size = 9),
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 10),
    axis.text.y = element_text(size = 10),
    axis.line = element_line(linewidth = 0.8),
    axis.ticks = element_line(linewidth = 0.8),
    panel.grid = element_blank(),        # removes all grid lines
    plot.title = element_text(size = 11, face = "bold", hjust = 0.5),
    plot.margin = margin(5, 5, 5, 5)
  )

# Save as high-quality PDF for publication
ggsave(
  "plot_interaction_meat_intentions.pdf",
  plot_interact_intentions_meat,
  height = 14, width = 20, units = "cm",
  dpi = 600, limitsize = FALSE
)

```


# Descriptive Statistics
```{r}
## Descriptive stats of DVs
summary(tot_df$quality_index)
sd(tot_df$quality_index)

summary(tot_df$health_index)
sd(tot_df$health_index)

summary(tot_df$convenience_index)
sd(tot_df$convenience_index)

summary(tot_df$environment_index)
sd(tot_df$environment_index)

summary(tot_df$animal_index)
sd(tot_df$animal_index)

summary(tot_df$price_index)
sd(tot_df$price_index)

table(tot_df$meat_intentions_1)
summary(tot_df$meat_intentions_1)
sd(tot_df$meat_intentions_1)

table(tot_df$meat_intentions_2)
summary(tot_df$meat_intentions_2)
sd(tot_df$meat_intentions_2)

table(tot_df$sub_intentions_1)
summary(tot_df$sub_intentions_1)
sd(tot_df$sub_intentions_1)

table(tot_df$sub_intentions_2)
summary(tot_df$sub_intentions_2)
sd(tot_df$sub_intentions_2)

table(tot_df$policy_support_1)
summary(tot_df$policy_support_1)
sd(tot_df$policy_support_1)

table(tot_df$policy_support_2)
summary(tot_df$policy_support_2)
sd(tot_df$policy_support_2)

table(tot_df$policy_support_3)
summary(tot_df$policy_support_3)
sd(tot_df$policy_support_3)

## Treatment and control variables
table(tot_df$group)
table(tot_df$group)/1219

summary(tot_df$age)
tot_df$age_group <- cut(
  tot_df$age,
  breaks = c(18, 25, 35, 45, 55, 65, Inf),
  labels = c("18-24", "25-34", "35-44", "45-54", "55-64", "65+"),
  right = FALSE
)
prop.table(table(tot_df$age_group)) * 100
table(tot_df$gender)/1219

table (tot_df$education)/1219

table (tot_df$diet)/1219

table(tot_df$left_right)
summary(tot_df$left_right)

table(tot_df$urban_rural_true)
summary(tot_df$urban_rural_true)

(table(tot_df$meat_consumption)/1219)*100

(table(tot_df$substitute_consum)/1219)*100

summary(tot_df$sub_general_opinion_w1)
sd(tot_df$sub_general_opinion_w1)

(table(tot_df$sub_brand_awareness_1_w1)/1219)*100

(table(tot_df$UserLanguage)/1219)*100
```



